import time

import numpy as np
import matplotlib

matplotlib.use('Agg')
from .mechine import Heat_recovery_system, steam_turbine, Condenser, D_pre, zero_pre
import threading
from sklearn.metrics import mean_squared_error, r2_score
import logging
from collections import deque
from pyswarm import pso
from logging.handlers import RotatingFileHandler
from influxDBturbinecooltest import Real_time_from_other_buffer_to_alg
import matplotlib.pyplot as plt
import os
import math

sim_steps = 0
ita_current = 0.85
kc_current = 1.38
# 7月参数
# stick_ita = 0.905
# stick_kc = 1.57
# 2月参数
# stick_ita = 0.72
# stick_kc = 2.42
# 3月参数
# stick_ita = 0.78
# stick_kc = 1.9
# ita_current = stick_ita
# kc_current = stick_kc
def share_data():
    global ita_current
    global kc_current
    global sim_steps
    return ita_current, kc_current, sim_steps

def setup_logger(thread_name):
    logger = logging.getLogger(thread_name)
    logger.setLevel(logging.DEBUG)  # 设置日志级别

    fmt = "[%(asctime)s-%(threadName)s] (%(filename)s %(lineno)d): %(levelname)s %(message)s"

    # 创建文件处理器
    log_filename = rf".\logs\{thread_name}.log"
    file_handler = RotatingFileHandler(log_filename, mode='a', maxBytes=5 * 1024 * 1024, backupCount=1, encoding="utf-8")
    file_handler.setLevel(logging.DEBUG)
    file_handler.setFormatter(logging.Formatter(fmt=fmt, datefmt="%Y-%m-%d %H:%M:%S"))
    logger.addHandler(file_handler)

    return logger

def optimize_ita(a, b, logger):
    """
    计算比例系数 k，使得 k * a 尽可能拟合 b
    :param a: 列表 a
    :param b: 列表 b
    :return: 比例系数 k
    """
    if len(a) != len(b):
        raise ValueError("列表 a 和 b 的长度必须相同")

    # 计算 k 的分子和分母
    numerator = sum(ai * bi for ai, bi in zip(a, b))
    denominator = sum(ai ** 2 for ai in a)

    # 避免除零错误
    if denominator == 0:
        raise ValueError("列表 a 的所有元素不能为零")

    # 计算比例系数 k
    k = numerator / denominator
    return k

# ========== 2. PSO 优化函数 for 凝汽器参数 ==========
def pso_optimize_kc(true_Pc_list, cond_input_list, kc_current, logger):
    """
    利用保存的数据辨识凝汽器模型中的流量（传热）系数参数。

    参数:
        true_Pc_list: [true_Pb, ...]，每个是真实背压（单位：KPa）
        cond_input_list: 每个元素为 [hc, Dc, Pb, n, Ta, h]，
                         其中 Pb 为传入 condenser.simulation 的背压输入
        kc_current: 当前凝汽器参数（初始值，默认为1.38）

    目标:
        调用 Condenser.simulation，用 candidate 参数替代原参数，
        得到预测背压；计算预测背压与真实背压之间的均方误差（MSE），最小化该误差。
    """

    def objective(x):
        candidate_param = x[0]
        sim_Pc = []
        condenser_model = Condenser()  # 创建新模型实例
        for vec in cond_input_list:
            hc, Dc, Pc_input, n, Ta, h = vec
            # 调用 condenser_model.simulation，使用 candidate_param 替代原参数
            # 返回值: (deta_Pc, Pc_sim, Dec, Tc, Hc1)
            _, Pc_sim, _, _, _, _, _, _ = condenser_model.simulation(hc, Dc, Pc_input, n, Ta, h, candidate_param)
            sim_Pc.append(Pc_sim)
        mse_val = mean_squared_error(true_Pc_list, sim_Pc)
        return mse_val

    lb = [0.75]
    #ub = [1.75]
    ub = [3]
    #ub = [10]
    try:
        xopt, fopt = pso(objective, lb, ub, swarmsize=50, maxiter=100, debug=False)
        return xopt[0]
    except Exception as e:
        logger.error(f"出现错误：{e}")
        return kc_current

def mechanics_core(logger):
    global ita_current
    global kc_current
    global sim_steps
    state_result = None
    max_len = 1500
    # simtime = 5000
    h = 0.1
    sample_time=5
    # sample_time=10
    tolerance=0.05
    max_delta_ita=0.05
    alpha=0.75
    data_save=0.01
    # data_save = 0.005
    max_delta_kc = 0.2
    # max_delta_kc = 1
    identify_start = 0.03
    # identify_start = 0.01
    show_progress=False

    turbine = steam_turbine(0.85)
    heatsystem = Heat_recovery_system()
    coldsystem = Condenser()

    Ne_sim = deque(maxlen=max_len)  # 用于计算误差值的仿真数据
    Ne_r = deque(maxlen=max_len)  # 用于计算误差值的真实数据
    Ne_DTs = []  # 用于辨识的真实数据
    Ne_sims = []
    turbine_inputs = []

    # 用于记录凝汽器的背压数据
    Pc_sim_list = deque(maxlen=max_len)  # 用于计算误差值的仿真背压（单位KPa）
    Pc_true_list = deque(maxlen=max_len)  # 用于计算误差值的真实背压（单位KPa）
    Cond_wat_D_true_list = deque(maxlen=max_len)
    Cond_wat_D_sim_list = deque(maxlen=max_len)
    Dec=0
    Pc_DTs = []  # 用于辨识的真实数据
    cond_inputs = []

    steady_steps = 0
    steady_state_reached = False

    ita_history = deque(maxlen=max_len)
    ita_history.append(ita_current)

    # 初始凝汽器参数
    kc_current_history = deque(maxlen=max_len)
    kc_current_history.append(kc_current)
    while True:
        # 从buffer中读取全部测点最新数据
        time.sleep(5)
        need_time, latest_data_other_new = Real_time_from_other_buffer_to_alg(logger)
        # input_points = [0, 1, 4, 2, 6, 7, 13, 12, 11, 3, 19, 14, 18]
        # data = latest_data_other_new[0, :][input_points]
        data = latest_data_other_new[0, :]
        logger.info(f'参数辨识输入数据{data}')
        data = data[1:]
        non_zeros = [0,1,2,3,4,5,6,7,10,11,12,13,14,17,18,19,22,23]
        if np.any(data[non_zeros] == 0):
            time.sleep(5)
            continue
        if state_result is None or len(state_result) == 0:
            state_result = [
                330, 3.571, 2.261, 2.1, 2, 480, 500, 520, 900, 1.203,
                0.6006, 0.3668, 0.1441, 0.04862, 62, 70, 77.85, 253, 5, 350,
                83, 95, 107, 5, 474, 113, 120, 138, 5, 712,
                167, 180, 188, 5, 823, 200, 210, 218, 5, 963,
                220, 236, 245, 54, 180, 0.01, 0.51, 0.00975
            ]
        # 解包当前状态
        (Ne, Ps1, Ps2, P1, P2, T1, T2, Tj, Ty2, Ps3, Ps4, Ps5, Ps6, Ps7,
         Tn7, Tl7, Tw7, hd7, Dd6, hd6, Tn6, Tl6, Tw6, Dd5, hd5, Tn5,
         Tl5, Tw5, Dd3, hd3, Tn3, Tl3, Tw3, Dd2, hd2, Tn2, Tl2, Tw2,
         Dd1, hd1, Tn1, Tl1, Tw1, Twc_m, Tw4, drhol4_dt, rhol4_m, Pc_sim
         ) = state_result

        # 读取当前时间步真实数据
        Ne_DT = data[0]  # 真实汽轮机负荷
        Pc_DT = data[17]  # 背压（单位：KPa）
        Pt_DT = data[3]
        Dsup_DT = data[1]
        Pre_real = data[5]
        Re_wat_T = data[6]
        Cond_wat_T = data[12]
        Cond_wat_D = data[11]
        Feed_wat_D = data[10]
        Sup_wat_T = data[2]
        n = data[13]
        Ta = data[18]

        Ne_r.append(Ne_DT)
        Ne_sim.append(Ne)

        # 保存真实背压（转换为MPa）
        Pc_true = Pc_DT / 1000
        Pc_true_list.append(Pc_DT)
        Cond_wat_D_true_list.append(Cond_wat_D)
        Pc_sim_list.append(Pc_sim * 1000)
        state_result_new = []
        # 内层仿真：调用汽轮机和凝汽器模拟
        for j in range(int(sample_time / h)):
            # 调用汽轮机模拟
            Ds, pt, ht, D1, p1, D1j, h1, D2, p2, Dj2, h2, P2, T2, D2r, \
                D3, p3, h3, Dj3, D4, p4, h4, Dj4, D5, p5, h5, Dj5, D6, p6, \
                h6, Dj6, D7, p7, h7, Dj7, hc, Dc, P1, T1, Tj, Ty2, Ne, \
                Cond_T, Ne8 = turbine.simulation(
                Sup_wat_T, Pt_DT, Dsup_DT, Ps1, Ps2, 0, Ne_DT, P1, P2,
                Pre_real, T1, T2, Re_wat_T, Tj, Ty2, Ps3, Ps4, Ps5, Ps6, Ps7,
                Pc_sim, h
            )

            # 使用当前 kc_current 更新凝汽器模拟
            deta_Pc, Pc_sim, Dec, Tc, Hc1, Vf, NTU, Qc = coldsystem.simulation(
                hc, Dc, Pc_sim, n, Ta, h, kc_current
            )
            # Pc_sim 为凝汽器模拟得到的背压（单位KPa）

            (Ps7, Dd7, hd7, Tn7, Tl7, Tw7, Ps6, Dd6, hd6, Tn6, Tl6, Tw6,
             Ps5, Dd5, hd5, Tn5, Tl5, Tw5, Ps4, Tw4, Ps3, Dd3, hd3, Tn3,
             Tl3, Tw3, Ps2, Dd2, hd2, Tn2, Tl2, Tw2, Ps1, Dd1, hd1, Tn1,
             Tl1, Tw1, rhol4, Twc_m, Tw4_m, rhol4_m, drhol4_dt) = heatsystem.simulation(
                Cond_wat_T, Tn7, Tl7, Ps7, Tw7, D7, h, hd7, Dd6, hd6,
                Tn6, Tl6, Ps6, Tw6, D6, h, Dd5, hd5, Tn5, Tl5, Ps5, Tw5, D5, h,
                Ps4, Dd3, hd3, D4, h, Dsup_DT, Tn3, Tl3, Ps3, Tw3, D3, h,
                Dd2, hd2, Tn2, Tl2, Ps2, Tw2, D2, h, Dd1, hd1, Tn1, Tl1, Ps1, Tw1,
                D1, h, Cond_wat_D, Feed_wat_D, Tw4, Twc_m, drhol4_dt, rhol4_m, h
            )

            state_result_new = [
                Ne, Ps1, Ps2, P1, P2, T1, T2, Tj, Ty2, Ps3, Ps4, Ps5, Ps6, Ps7,
                Tn7, Tl7, Tw7, hd7, Dd6, hd6, Tn6, Tl6, Tw6, Dd5, hd5, Tn5,
                Tl5, Tw5, Dd3, hd3, Tn3, Tl3, Tw3, Dd2, hd2, Tn2, Tl2, Tw2,
                Dd1, hd1, Tn1, Tl1, Tw1, Twc_m, Tw4, drhol4_dt, rhol4_m, Pc_sim
            ]
        if not (0 in state_result_new or any(math.isnan(x) for x in state_result_new)):
            state_result=state_result_new
        Cond_wat_D_sim_list.append(Dec)
        # Ne_sim.append(Ne)
        # Pc_sim_list.append(Pc_sim)  # 记录当前时间步的模拟背压

        # 计算汽轮机误差（以 Ne 为例）
        Ne_error = abs((Ne_DT - Ne) / Ne_DT)
        if Ne_error > data_save:
            turbine_input = [
                Sup_wat_T, Pt_DT, Dsup_DT, Ps1, Ps2, 0, Ne_DT, P1, P2,
                Pre_real, T1, T2, Re_wat_T, Tj, Ty2, Ps3, Ps4, Ps5, Ps6, Ps7,
                Pc_sim, h
            ]
            Ne_DTs.append(Ne_DT)
            Ne_sims.append(Ne)
            turbine_inputs.append(turbine_input)

        # 计算凝汽器误差（以背压为例），真实背压为 Pc_DT/1000
        Pc_error = abs((Pc_true - Pc_sim) / Pc_true)
        if Pc_error > 0.005:
            cond_input = [hc, Dc, Pc_sim, n, Ta, h]
            Pc_DTs.append(Pc_true)
            cond_inputs.append(cond_input)

        # 判定稳态（连续5次汽轮机误差小于 tolerance）
        if not steady_state_reached:
            if Ne_error < tolerance:
                steady_steps += 1
                if steady_steps >= 5:
                    steady_state_reached = True
                    logger.info("模型已进入稳态，开始参数辨识...")
            else:
                steady_steps = 0

        # 当汽轮机误差超过 identify_start 且已有样本，进行PSO辨识汽轮机效率
        if Ne_error > identify_start and len(turbine_inputs) > 0:
            new_ita = ita_current * optimize_ita(Ne_sims, Ne_DTs, logger)
            delta_ita = new_ita - ita_current
            delta_ita_limited = np.clip(delta_ita, -max_delta_ita, max_delta_ita)
            # ita_current = alpha * (ita_current + delta_ita_limited) + (1 - alpha) * ita_current
            ita_current = ita_current + alpha * delta_ita_limited
            #ita_current = stick_ita
            logger.info(f"汽轮机误差超过{identify_start * 100:.0f}%，PSO辨识效率：{ita_current:.4f}")
            turbine.ita = ita_current
            Ne_DTs.clear()
            Ne_sims.clear()
            turbine_inputs.clear()

        # 当凝汽器误差超过 tol_kc_identify 且已有样本，进行PSO辨识凝汽器参数
        if Pc_error > 0.02 and len(cond_inputs) > 0:
            kc_new = pso_optimize_kc(Pc_DTs, cond_inputs, kc_current, logger)
            delta_kc = kc_new - kc_current
            delta_kc_limited = np.clip(delta_kc, -max_delta_kc, max_delta_kc)
            # kc_current = alpha * (kc_current + delta_kc_limited) + (1 - alpha) * kc_current
            kc_current = kc_current + alpha * delta_kc_limited
            #kc_current = stick_kc
            logger.info(f"凝汽器误差超过{identify_start * 100:.0f}%，PSO辨识凝汽器参数：{kc_current:.4f}")
            Pc_DTs.clear()
            cond_inputs.clear()

        ita_history.append(ita_current)
        kc_current_history.append(kc_current)

        sim_steps = sim_steps+1

        # --- 仿真结束，输出并绘图 ---
        save_path = "pics/"
        os.makedirs(save_path, exist_ok=True)

        print(f"simsteps{sim_steps}")
        if sim_steps > 20:
        # if 20 < sim_steps < 100:
            print(f"simsteps{sim_steps}")
            Ne_sim_np = np.array(Ne_sim)
            Ne_r_np = np.array(Ne_r)
            mse = mean_squared_error(Ne_r_np, Ne_sim_np)
            MAPE = np.mean(np.abs((Ne_r_np - Ne_sim_np) / Ne_r_np)) * 100
            if len(Ne_r_np) > 1:
                R2 = r2_score(Ne_r_np, Ne_sim_np)
            else:
                R2 = float('nan')
            print(f"负荷评价指标：MSE={mse:.4f}, MAPE={MAPE:.4f}%, R2={R2:.4f}")

            # 绘制汽轮机效率变化曲线
            plt.figure(figsize=(10, 5))
            plt.plot(np.array(ita_history), color="b", label="汽轮机效率(估计值)")
            plt.xlabel("时间步", fontsize=10.5)
            plt.ylabel("效率", fontsize=10.5)
            plt.title("PSO辨识后汽轮机效率变化", fontsize=12)
            plt.legend()
            plt.savefig(f"{save_path}/汽轮机效率变化.png", dpi=300)
            plt.close()

            # 绘制凝汽器参数变化曲线
            if len(np.array(kc_current_history)) > 0:
                plt.figure(figsize=(10, 5))
                plt.plot(np.array(kc_current_history), marker='o', linestyle='-', color='purple', label="凝汽器流量系数")
                plt.xlabel("时间步", fontsize=10.5)
                plt.ylabel("流量系数", fontsize=10.5)
                plt.title("PSO辨识后凝汽器流量系数", fontsize=12)
                plt.legend()
                plt.savefig(f"{save_path}/凝汽器流量系数.png", dpi=300)
                plt.close()

            # 绘制负荷对比图
            plt.figure(figsize=(10, 5))
            plt.plot(Ne_sim_np, label="模拟负荷", color="g")
            plt.plot(Ne_r_np, label="真实负荷", color="r", linestyle='--')
            plt.xlabel("时间步", fontsize=10.5)
            plt.ylabel("负荷值", fontsize=10.5)
            plt.title(f"负荷对比\nMSE={mse:.4f}, MAPE={MAPE:.4f}%, R2={R2:.4f}", fontsize=12)
            plt.legend()
            plt.savefig(f"{save_path}/负荷比较.png", dpi=300)
            plt.close()
            print(f'模拟负荷值{Ne_sim_np}')

            # 绘制背压对比图
            Pc_true_array = np.array(Pc_true_list)
            Pc_sim_array = np.array(Pc_sim_list)
            mse_pc = mean_squared_error(Pc_true_array, Pc_sim_array)
            MAPE_pc = np.mean(np.abs((Pc_true_array - Pc_sim_array) / Pc_true_array)) * 100
            if len(Pc_true_array) > 1:
                R2_pc = r2_score(Pc_true_array, Pc_sim_array)
            else:
                R2_pc = float('nan')
            print(f"背压评价指标: MSE={mse_pc:.4f}, MAPE={MAPE_pc:.4f}%, R2={R2_pc:.4f}")

            plt.figure(figsize=(10, 5))
            plt.plot(Pc_sim_array, label="模拟背压", color="g")
            plt.plot(Pc_true_array, label="真实背压", color="r", linestyle='--')
            plt.xlabel("时间步", fontsize=10.5)
            plt.ylabel("背压 (KPa)", fontsize=10.5)
            plt.title(f"背压对比\nMSE={mse_pc:.4f}, MAPE={MAPE_pc:.4f}%, R2={R2_pc:.4f}", fontsize=12)
            plt.legend()
            plt.savefig(f"{save_path}/背压比较.png", dpi=300)
            plt.close()
            print(f'模拟背压值{Pc_sim_list}')

            # 绘制背压对比图
            CondD_true_array = np.array(Cond_wat_D_true_list)
            CondD_sim_array = np.array(Cond_wat_D_sim_list)
            plt.figure(figsize=(10, 5))
            plt.plot(CondD_sim_array, label="模拟凝结水流量", color="g")
            plt.plot(CondD_true_array, label="真实凝结水流量", color="r", linestyle='--')
            plt.xlabel("时间步", fontsize=10.5)
            plt.ylabel("凝结水流量", fontsize=10.5)
            plt.title(f"凝结水流量对比", fontsize=12)
            plt.legend()
            plt.savefig(f"{save_path}/凝结水流量比较.png", dpi=300)
            plt.close()
            #print(f'模拟凝结水流量值{CondD_sim_array}')

def param_ident():
    logger = setup_logger("param_ident")
    predict_thread = threading.Thread(target=mechanics_core, args=(logger,))
    predict_thread.start()